#include "csimplexutils.h"

/**
 * solveSimplex - finds optimal solution using two phase simplex method
 * @param a left sides of restrictions (matrix, input parameter)
 * @param b input right sides of restrictions (vector, input)
 * @param c input effectiveness function (vector, input)
 * @param type array of relations between left and right sides (input)
 * @param rc coefficients of resulting effectiveness function (vector, output)
 * @param rb coefficients of resulting right side (vector, output)
 * @param opt optimal value (output)
 * @return status
 */
template <typename T> int CSimplexUtils<T>::solveSimplex(const CMatrix<T> &a,
    const CVector<T> &b, const CVector<T> &c, int type[], CVector<T> &rc,
    CVector<T> &rb, T &opt) {

  // Left sides and right sides must have same number of rows
  if (a.getNumOfRows() != b.getSize()) {
    throw CInvalidSizeException(getClassName(), "solveSimplex");
  }

  // Left side and effectiveness function must have same number of columns
  if (a.getNumOfCols() != c.getSize()) {
    throw CInvalidSizeException(getClassName(), "solveSimplex");
  }

  // Hashmap which assigns simplex table column index to row index of ones in
  // basic columns
  map<TIndex, TIndex> identColIds;

  // Create simplex table and save index of column from which artificial
  // columns of table begins
  long firstFId(-1);
  CMatrix<T> table = createSimplexTable(a, b, c, type, identColIds, firstFId);

  // If simplex table contains artificial variables, perform first phase of
  // two phase simplex method (find initial acceptable solution) and pass
  // table to second phase
  if (firstFId != -1) {
    if (!doSimplex(table, identColIds, 2)) return SIMPLEX_SOLUTION_UNBOUNDED;
    prepareSecondSimplexPhase(table, firstFId);
  }

  // Second phase (find optimal solution)
  if (!doSimplex(table, identColIds, 1)) return SIMPLEX_SOLUTION_UNBOUNDED;

  // Find optimum, optimal effectiveness function and optimal right side
  // in resulting simplex table
  return getSolutionFromTable(table, identColIds, rc, rb, opt);
}
//---------------------------------------------------------------------------

/**
 * createSimplexTable - creates simplex table
 * @param a left sides of restrictions (matrix, input)
 * @param b right sides of restrictions (vector, input)
 * @param c effectiveness function (vector, input)
 * @param type array of relations between left and right sides (input)
 * @param identColIds indices of identity columns (input)
 * @param index of first artificial column (output)
 * @return matrix representing simplex table (output)
 */
template <typename T> CMatrix<T> CSimplexUtils<T>::createSimplexTable(
    const CMatrix<T> &a, const CVector<T> &b, const CVector<T> &c, int type[],
    map<TIndex, TIndex> identColIds, long &firstFId) {

  // Matrix for resulting simplex table and temporary vector for right side
  CMatrix<T> table(a);
  CVector<T> right(b);

  // 1) Adding auxiliary variables: slack, surplus and artificial
  //-------------------------------------------------------------------------

  // Temporary containers "sCols" for slack and surplus variables, "aCols"
  // for artificial variables
  vector<CVector<T> > sCols, aCols;
  TIndex auxColCount = 0;

  // Analyze restrictions
  for (TIndex i = 0; i < table.getNumOfRows(); i++) {

    // If right side of actual restriction is negative, multiply both sides by
    // -1
    if (right.get(i) < 0) {
      right.set(i, -right.get(i));
      for (TIndex j = 0; j < table.getNumOfCols(); j++) {
        table.set(i, j, -table.get(i, j));
      }
      type[i] *= -1;
    }

    // If actual restriction contains "<=" add surplus variable to "sCols"
    // If actual restriction contains "=" add artificial variable to "aCols"
    // If actual restriction contains ">=" add slack (negative) variable
    // "sCols" and artificial variable to "aCols"
    switch (type[i]) {
      case CONDITION_LEFT_SIDE_LESSER: {
        sCols.push_back(createAuxCol(table.getNumOfRows(), i, true));
        identColIds[i] = table.getNumOfCols() + auxColCount++;
      } break;
      case CONDITION_LEFT_SIDE_EQUAL: {
        aCols.push_back(createAuxCol(table.getNumOfRows(), i, true));
        identColIds[i] = table.getNumOfCols() + auxColCount++;
      } break;
      case CONDITION_LEFT_SIDE_GREATER: {
        sCols.push_back(createAuxCol(table.getNumOfRows(), i, false));
        aCols.push_back(createAuxCol(table.getNumOfRows(), i, true));
        identColIds[i] = table.getNumOfCols() + (++auxColCount)++;
      } break;
      default: break;
    }
  }

  // Add surplus/slack and artificial columns to the table (if exist)
  if (!sCols.empty()) {
    table = CMatrix<T>(table, CMatrix<T>(CMatrix<T>::VECTOR_COL, sCols));
  }
  if (!aCols.empty()) {
    firstFId = (long) table.getNumOfCols();
    table = CMatrix<T>(table, CMatrix<T>(CMatrix<T>::VECTOR_COL, aCols));
  }

  // 2) Adding auxiliary effectiveness function
  //-------------------------------------------------------------------------

  // Add right side and effectiveness function to the table
  table = CMatrix<T>(table, right, CMatrix<T>::VECTOR_COL);
  table = CMatrix<T>(table, c, CMatrix<T>::VECTOR_ROW);

  // If table contains artificial variables, create auxiliary effectiveness
  // function
  if (firstFId != -1) {

    // Auxiliary effectiveness function for first phase
    CVector<T> tmpC(table.getNumOfCols());

    // Iterate through columns of simplex table
    for (TIndex i = 0; i < table.getNumOfCols(); i++) {

      // Skip artificial column
      if (i >= (TIndex) firstFId && i < table.getNumOfCols() - 1) continue;

      // Temporary sum of items in actual column
      T sum(0);

      // Create sum from those items in actual column, that are located in
      // rows, that contain number one in one of artificial columns
      // Actual item of auxiliary effectiveness function (under current
      // column) will contain opposite value of this sum
      for (TIndex j = 0; j < table.getNumOfRows(); j++) {
        bool hasArtificial = false;
        for (TIndex k = (TIndex) firstFId; k < table.getNumOfCols() -1; k++) {
          if (table.get(j, k) == 1) {
            hasArtificial = true;
            break;
          }
        }
        if (hasArtificial) sum += table.get(j, i);
      }
      tmpC.set(i, -sum);
    }

    // Add auxiliary effectiveness function to the table
    table = CMatrix<T>(table, tmpC, CMatrix<T>::VECTOR_ROW);
  }
  return table;
}
//---------------------------------------------------------------------------

/**
 * createAuxCol - create auxiliary (surplus/slack/artificial)
 * @param size number of items in column
 * @param onePosition position of one in identity column
 * @param onePositive true if number one will be positive
 * @return auxiliary column
 */
template <typename T> CVector<T> CSimplexUtils<T>::createAuxCol(TIndex size,
    TIndex onePosition, bool onePositive) {
  CVector<T> tmp(size);
  tmp.clearToZero();
  tmp.set(onePosition, onePositive ? 1 : -1);
  return tmp;
}
//---------------------------------------------------------------------------

/**
 * doSimplex - performs simplex algorithm
 * @param table simplex table
 * @param identColIds indices of identity columns
 * @param cRowCount number of effectiveness functions
 * @return true if optimal solution has been found
 */
template <typename T> bool CSimplexUtils<T>::doSimplex(CMatrix<T> &table,
    map<TIndex, TIndex> identColIds, TIndex cRowCount) {

  // While optimal solution is not found or problem is unbounded, do
  // elimination
  while (true) {

    // Row and column of pivot
    TIndex pivotRow, pivotCol;

    // Optimality test, elimination
    switch (classifyTable(table, identColIds, cRowCount, pivotRow, pivotCol)) {
      case BLAND_SOLUTION_NONOPTIMAL: table.elimCol(pivotRow, pivotCol); break;
      case BLAND_SOLUTION_OPTIMAL: return true;
      case BLAND_SOLUTION_UNBOUNDED: return false;
      default: break;
    }

    // Change index of identity (base) column
    identColIds[pivotRow] = pivotCol;
  }
}
//---------------------------------------------------------------------------

/**
 * classifyTable - performs optimality test, Bland's rule
 * @param table simplex table
 * @param cRowCount number of effectiveness functions
 * @param pivotRow index of pivot's row (output)
 * @param pivotCol index of pivot's column (output)
 * @return result of Bland's rule (solution, opt. solution, unbounded)
 */
template <typename T> int CSimplexUtils<T>::classifyTable(CMatrix<T> &table,
    map<TIndex, TIndex> identColIds, TIndex cRowCount, TIndex& pivotRow,
    TIndex& pivotCol) {

  // Temporary indices
  TIndex lastRowId(table.getNumOfRows() - 1);
  TIndex lastColId(table.getNumOfCols() - 1);
  TIndex lastButOneColId(table.getNumOfCols() - 2);
  TIndex numOfLeftSideCols(table.getNumOfCols() - 1);

  // Set pivot's row and column index to zero
  pivotRow = 0;
  pivotCol = 0;

  // Bland's rule
  //-------------------------------------------------------------------------
  // pivotCol = s = min{j; c[j] < 0}
  // pivotRow = r pro ktere B[r] = min{B[k]; b[k] / a[k][s] = min{
  // b[j] / a[j][s]; a[j][s] > 0}, a[k][s] > 0}
  //-------------------------------------------------------------------------

  // Get index of first negative item in efectiveness function (except right
  // side), if there is no such item, table is optimal, return
  for (TIndex i = 0; i < numOfLeftSideCols; i++) {
    T tmp(table.get(lastRowId, i));
    if (tmp < 0) {
      pivotCol = i;
      break;
    }
    if (i == lastButOneColId) return BLAND_SOLUTION_OPTIMAL;
  }

  // Get subtable, which does not contain efectiveness functions
  CMatrix<T> tmpTable(table.subMatrix(0, 0, table.getNumOfRows() - cRowCount,
      table.getNumOfCols()));

  // Temporary zero vector
  CVector<T> zero(tmpTable.getNumOfRows());
  zero.clearToZero();

  // If subtable has non-positive column above first negative item in
  // efectiveness function, problem is unbounded, return
  if (tmpTable.getCol(pivotCol) <= zero) return BLAND_SOLUTION_UNBOUNDED;

  // Find row in which is right side item (in last column) divided by item
  // in pivot's column minimal
  T tmpDiv(0);
  bool firstMatch(true);
  for (TIndex i = 0; i < tmpTable.getNumOfRows(); i++) {

    // Calculate temporary quotient
    if (tmpTable.get(i, pivotCol) <= 0) continue;
    T tmp = tmpTable.get(i, lastColId) / tmpTable.get(i, pivotCol);

    // If actual temporary quotient is first quotient or is smallest from all
    // found quotients, save it's index
    // In case that two smallest equal quotients are found, take that in
    // which row is number one in most left identity (base) column
    if (tmp < tmpDiv || firstMatch) {
      pivotRow = i;
      tmpDiv = tmp;
      firstMatch = false;
    } else if (tmp == tmpDiv && identColIds[i] < identColIds[pivotRow]) {
      pivotRow = i;
      tmpDiv = tmp;
    }
  }
  return BLAND_SOLUTION_NONOPTIMAL;
}
//---------------------------------------------------------------------------

/**
 * prepareSecondSimplexPhase - removes auxiliary function and columns
 * @param table table
 * @param firstFId index of first artificial column
 */
template <typename T> void CSimplexUtils<T>::prepareSecondSimplexPhase(
    CMatrix<T> &table, long firstFId) {

  // Indices of last row and column
  TIndex lastRowId = table.getNumOfRows() -1;
  TIndex lastColId = table.getNumOfCols() -1;

  // Remove auxiliary efectiveness function
  table = table.subMatrix(0, 0, lastRowId, table.getNumOfCols());

  // Remove artificial columns (they are starting at "firstAId", but behind
  // them is right side column)
  CVector<T> tmpRight = table.getCol(lastColId);
  table = table.subMatrix(0, 0, table.getNumOfRows(), firstFId);
  table = CMatrix<T>(table, tmpRight, CMatrix<T>::VECTOR_COL);
}
//---------------------------------------------------------------------------

/**
 * getSolutionFromTable - gets solution from resulting simplex table
 * @param table
 * @param identColIds
 * @param rc
 * @param rb
 * @param opt
 * @return information about number of solutions
 */
template <typename T> int CSimplexUtils<T>::getSolutionFromTable(
    CMatrix<T> &table, map<TIndex, TIndex> identColIds, CVector<T> &rc,
    CVector<T> &rb, T &opt) {

  // Indices of last row and column
  TIndex lastRId = table.getNumOfRows() -1;
  TIndex lastCId = table.getNumOfCols() -1;

  // Gets optimum
  opt = -table.getCol(lastCId).get(lastRId);

  // Gets resulting efectiv. function's row "rc" and right side column "rb"
  rc = table.subMatrix(0, 0, table.getNumOfRows(), lastCId).getRow(lastRId);
  rb = table.subMatrix(0, 0, lastRId, table.getNumOfCols()).getCol(lastCId);

  // If there is (rc[i] > 0 for every i that is not in identity column) in
  // resulting vector, then exists only one optimal solution
  for (TIndex i = 0; i < rc.getSize(); i++) {
    map<TIndex, TIndex>::iterator iter = identColIds.begin();
    for (; iter != identColIds.end(); ++iter) {
      if (i != iter->second && rc.get(i) == 0) return SIMPLEX_SOLUTION_INF;
    }
  }
  return SIMPLEX_SOLUTION_ONE;
}
//---------------------------------------------------------------------------

template class CSimplexUtils<char>;
template class CSimplexUtils<short>;
template class CSimplexUtils<int>;
template class CSimplexUtils<long>;
template class CSimplexUtils<float>;
template class CSimplexUtils<double>;
template class CSimplexUtils<CBigInt>;
template class CSimplexUtils<CFraction>;
